##Anderson-Nilsson & Clayton 
##Figures 1 & 2

rm(list=ls())
library(car)
library(psych)
library(texreg)
library(foreign)
library(readstata13)
library(ggplot2)
library(plyr)

# set your working directory here


data<-read.dta13('Persuasion_Study.dta')

##Figure 1
##just birth control 
data$agreement_bc_susan<-recode(data$agreement_bc_brian, "1=0; 0=1")

data4<-subset(data, treat7 ==1 | treat8 ==1)

data4$diff<-as.numeric(data4$susan_expertise_birthcontrol) - as.numeric(data4$brian_expertise_birthcontrol)


## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {

    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

    # This does the summary. For each group's data frame, return a vector with
    # N, mean, and sd
    datac <- ddply(data, groupvars, .drop=.drop,
      .fun = function(xx, col) {
        c(N    = length2(xx[[col]], na.rm=na.rm),
          mean = mean   (xx[[col]], na.rm=na.rm),
          sd   = sd     (xx[[col]], na.rm=na.rm)
        )
      },
      measurevar
    )

    # Rename the "mean" column    
    datac <- rename(datac, c("mean" = measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}

data_sum <- summarySE(data4, measurevar="agreement_bc_susan", groupvars=c("diff","party_id"))

data_sum<-subset(data_sum, party_id=="Democrat" | party_id == "Republican")


# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.1) # move them .05 to the left and right

breaks<-seq(-3,3,1)

pdf(file = "Figure1.pdf", 
    width = 6, 
    height = 6) 

ggplot(data_sum, aes(x= diff, y= agreement_bc_susan, colour= party_id, group= party_id)) + 
    geom_errorbar(aes(ymin= agreement_bc_susan-se, ymax= agreement_bc_susan +se), width=.1, position=pd) +
    geom_line(position=pd) +
    geom_point(position=pd, shape=21, fill="grey", aes(size=N), show.legend=FALSE) + # 21 is filled circle
    xlab("Difference in Expertise Evaluations (Susan - Brian)") +
    ylab("Percent Agree with Susan") +
    scale_colour_grey(name="Respondent Party",    # Legend label, use darker colors
                     breaks=c("Republican", "Democrat"),
                     labels=c("Republican", "Democrat"),
                     start = 0, end = 0.7,
    			    na.value = "grey50") +
    ggtitle("Average Agreement with Susan by Difference in Expertise Rating") +
    expand_limits(y=0:1) +                        # Expand y range
    scale_x_continuous(breaks= breaks ) +         # Set tick every 4
    theme_bw() +
    theme(legend.position=c(0.1,0.85)) +               # Position legend in bottom right
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
   panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(hjust=0.5))+
   geom_hline(yintercept=0.5, color="grey") +
  geom_vline(xintercept=0, color='grey') 

dev.off()

#Figure 2

##Stacked means 

#create four types: D/R m & w 

data4$GP<-paste(data4 $gender1, data4 $party_id, sep=" ")
data4$GP<-factor(data4 $GP, levels= c("Female Democrat", "Male Democrat", "Female Republican", "Male Republican"))

data4$BrianGen<-paste(data4 $gender1, data4 $BrianFirstvSusanFirst, sep=" ")
data4$BrianGen<-as.factor(data4 $BrianGen)


data_sum2 <- summarySE(data4, measurevar="favor_oppose_otcbc", groupvars=c("BrianFirstvSusanFirst", "party_id", "gender1"))
data_sum3$BrianFirstvSusanFirst <- factor(data_sum3$BrianFirstvSusanFirst)

data_sum2<-subset(data_sum2, party_id=="Democrat" | party_id=="Republican")

data_sum3 <- summarySE(data4, measurevar="favor_oppose_otcbc", groupvars=c("BrianFirstvSusanFirst", "GP"))

data_sum3<-subset(data_sum3, GP == "Female Democrat" | GP == "Male Democrat" | GP == "Female Republican" | GP == "Male Republican")

data_sum3$GP<-recode(data_sum3$GP, "'Female Democrat'= 'Democrats Women'; 'Male Democrat' = 'Democrats Men'; 'Female Republican'='Republicans Women'; 'Male Republican'='Republicans Men'")

data_sum3$BrianFirstvSusanFirst<-recode(data_sum3$BrianFirstvSusanFirst, "0='Susan'; 1='Brian'")

data_sum3$BrianFirstvSusanFirst<-as.factor(data_sum3$BrianFirstvSusanFirst)

levels(data_sum3 $GP) <- gsub(" ", "\n", levels(data_sum3 $G))


pdf(file = "Figure2.pdf", 
    width = 6, 
    height = 6) 

ggplot(data_sum3, aes(x= GP, y= favor_oppose_otcbc, fill= BrianFirstvSusanFirst)) + 
    geom_bar(position=position_dodge(), stat="identity") +
    geom_errorbar(aes(ymin= favor_oppose_otcbc-se, ymax= favor_oppose_otcbc +se),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+
    theme_bw() +
    xlab("") +
    ylab("Percent Supporting OTC Birth Control") +
    ggtitle("Average Agreement with OTC Birth Control by Speaker in Favor") +
    theme_bw() +    
    expand_limits(y=0:1) +                        # Expand y range
    theme(legend.position=c(0.87,0.9)) +               # Position legend in bottom right
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
   panel.background = element_blank(), axis.line = element_line(colour = "black"), plot.title = element_text(hjust=0.5))+    
    #scale_fill_discrete(name = "Speaker", labels = c("Susan", "Brian"))+
       scale_fill_manual("Speaker in Favor", values = c("Susan" = "grey40", "Brian" = "grey80")) +
       theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))

# Error bars represent standard error of the mean

dev.off()



